Stochastic Modelling Approach to the Incubation Time of Prionic Diseases 



A. S. Ferreirai, M. A. da SilvaB and J. C. Cressoni^l 

^ Universidade Federal de Alagoas, Departamento de Fisica, 57072-970 Maceio (AL), Brazil 
* Departamento de Fisica e Qmmica, FCFRP, Universidade de Sao Paulo, 14040-903 Riheirdo Preto, SP, Brazil 

(Dated: February 6, 2008) 

Transmissible spongiform encephalopathies like the bovine spongiform encephalopathy (BSE) and 
the Creutzfeldt- Jakob disease (CJD) in humans are neurodegenerative diseases for which prions are 
the attributed pathogenic agents. A widely accepted theory assumes that prion replication is due to 
a direct interaction between the pathologic (PrP^'^) form and the host encoded (PrP*^) conformation, 
in a kind of an autocatalytic process. Here we show that the overall features of the incubation time 
of prion diseases are readily obtained if the prion reaction is described by a simple mean-field model. 
An analytical expression for the incubation time distribution then follows by associating the rate 
constant to a stochastic variable log normally distributed. The incubation time distribution is then 
also shown to be log normal and fits the observed BSE data very well. The basic ideas of the 
theoretical model are then incorporated in a cellular automata model. The computer simulation 
results yield the correct BSE incubation time distribution at low densities of the host encoded 
protein. 

PACS numbers; 87.10.+e 87.19.Xx 05.20.Dd 



X 



The so-called prion diseases comprise a group of 
fatal transmissible spongiform encephalopathies (TSE) 
like the well known bovine spongiform encephalopathy 
(BSE) and sheep scrapie. In humans, these progressive 
neurodegenerative diseases include Kuru, Creutzfeldt- 
Jakob disease (CJD), Gerstmann-Straeussler-Scheinker 
syndrome (GSS) and fatal familial insomnia (FFI). Com- 
mon pathology includes spongiform degeneration and 
characteristic formation of plaques in the brain tissue Q . 
Variant CJD correlated with a (BSE)-like prion strain, 
have been identified and are believed to be linked to the 
consumption of contaminated food IE Q Hi - 

The protein-only hypothesis states that the infec- 
tious agent is a protein, named prion 0, ll| , which is de- 
void of nucleic acid and capable of replicating itself in the 
absence of these traditional genetic material. Two con- 
formations of this protein are important for characteriz- 
ing the disease, namely, the normally folded host-encoded 
cellular protein called PrP*^ and an abnormal pathogenic 
conformation named PrP^'^. The latter form is hydropho- 
bic, has a tendency to form aggregates and may be found 
in different strains. The pathogenic form PrP^'^ is more 
stable than the endogenous cellular form, and is known 
to be resistant to enzymatic digestion, radiation and high 
temperatures. One of the most accepted models for prion 
replication assumes that this form acts as a template for 
converting the host prion into its own conformation in a 
kind of an autocatalytic reaction 0| . Understanding 
the dynamics of the PrP*^ PrP^'^ transformation is 
crucial if one is attempting to explain and predict the 
main stages of the disease. The reaction is complex, 
perhaps involving other participants possibly acting as 
chaperone, to help mediate protein folding The 
number of parameters involved for thoroughly describing 
the transformation process is thus expected to be very 
large 0, 0| . It is therefore important to be able to rec- 



ognize which ones are mandatory, i.e., responsible for the 
major aspects of the dynamics. 

Here we present a simple, analytically solvable, mean- 
field model for describing the prion reaction problem, 
which focuses on realistically reproducing the incubation 
time of the disease. For notational convenience it is use- 
ful to introduce the following definitions: A stands for 
the host protein (PrP*~") and B stands for the pathogenic 
form (PrP'^'') with a — [A] and b — [B] denoting volume 
concentrations. We then write the autocatalytic conver- 
sion reaction simply as 



A + B 
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where K is the reaction rate. For simplicity we shall as- 
sume that the total concentration a + b ^ p is kept fixed 
at all times. This means that there is no metabolic de- 
composition of B and any metabolic decomposition of A 
is immediately compensated by the host genetic system. 
It also implies that the host takes no action for produc- 
ing new, normal protein, as the reaction takes place. In 
order to stick to the simplest possible case we are also as- 
suming that the reaction is unidirectional and favors the 
most stable form PrP^'^. No other strains are supposed to 
be present and both forms are assumed to be uniformly 
distributed. The kinetic evolution is then given by 
db/dt = Kab — K{p — b)b which is the simplest possible 
nonlinear equation describing an autocatalytic reaction. 
This equation can be easily integrated up to the time T 
giving 



T = 
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with bo being the infection dose given at time t — and 
ao the initial concentration of A. According to this ex- 
pression b(t) is slowly varying for small t, followed by 
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a period of rapid increase in a short time interval, then 
reaching a plateau for long enough times when the reac- 
tion stops [Hfisj. 

We now define the incubation time (Tj) as the time 
it takes for the number of pathogenic prions to reach a 
given value 6/, i.e., b{Ti) — hi. (It makes no difference to 
our calculations whether 6/ represents a number of prions 
or an aggregate with size 6/.) A useful approximation can 
be obtained by assuming, reasonably, that 6o/ao << 1. 
This gives 
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This log-dependence of the incubation time on the ini- 
tial dose was quantitatively observed by Prusiner jl^ 
from inoculation of a form of scrapie in hamsters (Fig.[T]l. 
Prusiner 's results also indicate that the survival time is 
practically independent of the dose. Eqn. Q is consis- 
tent with this finding (see also If we define the 
time of death as the time it takes for the number of 
B' s to reach the value h^, i.e., h{T£,) — bjj, we find 
that Tjj — Tj does not depend on bo- Moreover, eqn. © 
can be easily adapted to fit Prusiner's data. In order to 
mimic the end-point titration method used in the exper- 
iment, we first define all concentrations relative to the 
largest experimental concentration which we shall call 
Po- We then write bo/Po = 10"^^" (n = dose) and al- 
low n to vary from n = (smallest concentration) to 
n = 10 (largest concentration). We can now apply re- 
gression to the data (using only the integral values for 
n) to obtain the best fit. Notice, however, that the ex- 
perimental curves are composed of two branches, both 
exhibiting a sudden increase in the inclination for n < 2 
(see Fig. QJ. This behavior seems to be indicative of 
a threshold, possibly leading to a smaller rate constant 
at high dilutions. One can simulate a dose dependent 
activation mechanism linked to the rate constant K by 
making the following "ansatz" : we make K — > K^g with 
Keff ^ K[l — ai/{a2 + exp(n))]. With these implemen- 
tations, eqn. Q reads Ti = C — [In 10/(i4rao)]n, C be- 
ing a constant (independent of ho). The phenomenolog- 
ical constants, estimated with a non-linear least-squares 
fitting to this equation, with K replaced by Keff , are 
ai = 0.23(4) (0.61(2)) and aa = -0.51(6) (2.1(2)) for the 
incubation (death) curve. The result of the full fitting 
is shown in the inset of Fig. ^ Notice that Keff rapidly 
approaches K for n > 2. 

However, we decided to avoid dealing with the contro- 
versial features associated with the region n < 2 (contain- 
ing only two points) and stick to the (larger) less inclined 
part of the experimental curve. Therefore any param- 
eter obtained from the ?/— intercepts in Fig. ^ will not 
be taken into account. The regression coefficient gives 
l/{KaQ) = 3.12(3) days for the incubation part of the 
curve and \/{KaQ) = 3.02(6) days for the death part 



of the curve. This (partial) fitting is represented by the 
continuous line in the main part of Fig. We can easily 
check the reasonableness of these figures. Notice that we 
could have started with the Michaelis-Menten equation, 
namely, dh/dt = Kxlab/ {Km + a)] with Kt and Km be- 
ing the turnover number and the Michaelis constant re- 
spectively |l3j . Direct integration of this equation yields 
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which is consistent with eqn. (jSJ for Km » ao » bo 
and K ~ Kt/Km- We can therefore estimate the prion 
Kt/Km ratio for the scrapie strain used by Prusiner in 
hamsters. If we assume ao nanomole liter"^ |l2l fl^ 
we find Kt/Km ~ 10^ M~^s~^. This value is within 
the range expected for enzymes, in which case Km lies 
between 10~^ M to 10""'^ M and Kt falls in the range 
from 10 s^^ and 10^ s^^. 

Having discussed the behavior of the incubation time 
on bo, we now turn our attention to the dependence of 
T/ on ao- The role played by the host prion initial con- 
centration is useful for describing reactions, such as Ql, 
in numerical simulation approaches. The explicit power 
law dependence of T/ on ao can be seen by expanding 
eqn. |(5J in terms of bj/ao- This gives 
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and therefore T/ ^ Ai/ ao^ A2/ Oq- This kind of behavior, 
having the form of a sum of monomer and dimer terms, 
has already been suggested in the literature How- 
ever, the determination of the explicit dependence of the 
coefficients Ai and A2 on 6/ and bo, as shown here, was 
only possible because of the simplicity of the model. 

The initial concentration of the endogenous PrP pro- 
tein is determinant for the dynamics of the prion reaction 
since it represents the reaction fuel. The higher the initial 
concentration ao, the lower the time for the PrP^'^ con- 
centration to reach the value 6/ . These results have been 
obtained through careful computer simulations by Cox 
et al. They also showed that the incubation time 

distributions for different oq collapse to a single form 
if the time scale is properly normalized to unity. Will 
our simple, minimally parametrized model represented 
by the basic reaction QJ, be able to reproduce such re- 
sults? In order to address this question, we ran computer 
simulations based on a cellular automata {CA) with rules 
following a close resemblance to our model. 

According to the CA rules, a,n NxN {N — 200) square 
lattice is randomly populated with a number (Nao) of the 
host A = PrP^ protein and a number {Nbo = 6) of the 
B = PrP^'^ misfolded protein. Nao is given as a small 
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FIG. 1: Dependence of the incubation time (Tj) on the in- 
fection (initial) dose (= n with 60 //3o = 10""^°). The ex- 
perimental data were obtained from Prusiner's work and 
dashed lines are just meant to lead the eye. In the main figure 
we apply regression to the data (n > 2) to obtain the best fit 
with eqn. I^J . The most diluted part {n < 2) was left out due 
to the abrupt change in behavior in this region, leaving only 
two points (n = 1,2) for the fitting. Therefore only the incli- 
nations (= l/{Kao)) are kept. The inset shows a non-linear 
least-squares full fitting (all n) to eqn. ^ with the ansatz 
K Keff = K[l - ai/{a2 + exp(n))]. 



percentage of the total number of sites available and to 
each of the B sites is assigned a "mass" (m) , initially set 
to unity. The A's and B's are allowed to diffuse randomly 
to their nearest neighbor sites and a reaction occurs when 
a -B is approached by an A at a distance d < \pm. In this 
case the normal Prion disappears and the misfolded Prion 
has its mass increased by one. The reaction is unidirec- 
tional, favoring i?, with the A's slowly disappearing from 
the system, keeping A + B ^ constant. One site- by-site 
sweep through the lattice is made for diffusion followed by 
another one for reaction. The time unit is then increased 
by one (arbitrary units). The reaction stops when one of 
the masses reaches the value m = 40, the corresponding 
computer time thus characterizing the incubation time. 
The above values for the parameters (not the CA rules) 
were adjusted from the numerical simulations of Cox et 
al. [T^ in a hexagonal lattice. The mass is here to mimic 
clusterization without assigning any geometric form to 
the cluster. Besides simplifying the computer code and 
speeding up the simulations, this helps reducing the in- 
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FIG. 2: Incubation time distributions with the time scale 
normalized by the mean incubation time, i.e., tj — Ti/Tj. 
The full circles represent the observed incubation time distri- 
bution for BSE- infected cattle in UK |l3, [ll [l^ . Fig. Ha) 
shows the results from computer simulations on an A*' x A*' 
(A'' = 200) square lattice, with a number of PrP^'^ seeds 
Nbo = 6 (see text). The values of Nao shown represent 
the initial PrP*^ concentration and only a few curves for A^ao 
were drawn to avoid figure cluttering. Notice the tendency 
for better agreement with the observed results as Nao gets 
smaller. Fig. |5Jb) shows the same experimental data as in 
(a), along with the proposed analytical distribution G{ti), 
obtained from our model assuming a log normal distribution 
for the rate constant. When the time units are scaled by the 
mean time we are left with a single parameter, namely a, 
whose best fitted value is given by cr = 0.255 ± 6. 



fluence of local topology on the final results. 

Fig. |2Ia) shows the simulation results for the incuba- 
tion time distributions for several values of A'^o, with the 
time scale normalized by the mean time. Notice that as 
the PrP^ concentration is decreased, the corresponding 
distribution converge asymptotically to the experimen- 
tal results (BSE-infected cattle in UK [Ulilllil) rep- 
resented by the full circles. Increasing A^^o makes the 
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system more homogeneous which diminishes fluctuations 
and narrows the distribution. The biological concentra- 
tions (believed to be nanomolar) correspond to an areal 
concentration around N'^g = 0.001% [13 • With the CA 
rules above, such small concentrations would require very 
large computing time, if feasible at all. The best agree- 
ment is obtained for Nao — 0.6% which is far as we could 
go with these simulations. 

Our next issue is to search for an analytical form for the 
incubation time distribution. Knowledge of such a func- 
tion is not only important to check the reliability of the 
model but also to provide a distribution that can be used 
in statistical studies . We need to adapt the determin- 
istic model to accommodate a stochastic variable follow- 
ing a known distribution and associate it with eqn ^ 
(or Q). Since the protein-folding process actually in- 
volve many steps 0, possibly chaperone-assisted [2]| . 
the end result of the prionic reaction can adequately be 
viewed as a series of multiplicative processes. It is there- 
fore reasonable to assume that the distribution of the 
reaction rate K in a population is log normal |22j |. Since 
oc l/Tf it is easy to show that Tj also follows a log 
normal distribution with the same deviation. The scaled 
distribution G{ti), with tj ~ Tj/Tj, is then readily ob- 
tained. One finds 
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Gitj) = 



1 



tj exp 



Ini/ + {5/2)a' 



which does not depend neither on the initial variables 
nor on bj. We are therefore left with a single fitting pa- 
rameter, namely a, the standard deviation of InK. Ap- 
plying non-linear least squares fitting to Eqn. © we get 
a = 0.255(±6). The final result is shown in Fig.EJb). In 
this Figure the observed data are the same used in refer- 
ence 17] for BSE- infected cattle in the United Kingdom 
born in 1987 

In conclusion, a simple mean field model, based on an 
autocatalytic mechanism, is shown to contain the basic 
ingredients necessary to describe the essential features 
associated with the incubation time of the complex prion 
conversion reactions. Assuming that the rate constant 
is a random variable, following a log normal distribu- 
tion, we were able to provide a closed form for the in- 
cubation time distribution of BSE-infected cattle. The 
surprisingly simple analytical expression derived for the 
incubation time distribution contains only one param- 
eter, namely the variance of the logarithm of the rate 
constant. The simplicity of the model is characterized 
by the almost naive differential equation upon which it is 
based, by simple computer simulations, and by the mini- 
mal set of parameters used to describe the autocatalytic 
process. 
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